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Abstract 

The SOR iteration for solving linear systems of equations depends upon an overrelax- 
ation factor oj. We show that for the standard model problem of Poisson’s equation on a 
rectangle, the optimal u) and corresponding convergence rate can be rigorously obtained by 
Fourier analysis. The trick is to tilt the space-time grid so that the SOR stencil becomes 
symmetrical. The tilted grid also gives insight into the relation between convergence rates 
of several variants. 
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1. Introduction 


Fourier analysis has been used for nearly fifty years to test the stability of time- 
dependent finite difference formulas — the “von Neumann method” [9]. More recently it 
has also become a standard tool for estimating the convergence rate of multigrid iterations 
[3], But the analysis of the more classical iteration known as successive overrelaxation — 
SOR — has been carried out by other means [5,6,10,11]. The reason is that the behavior 
of SOR, unlike the other two problems, is dominated by low-frequency modes that are 
controlled by boundary conditions. The obvious application of Fourier analysis treats 
the boundary conditions incorrectly, and leads to an incorrect prediction of the optimal 
convergence rate. 

In this note we show that if the computational grid is tilted by a certain angle in 
space and time, Fourier analysis becomes exact for the standard SOR model problem: 
the five-point discretization of Poisson’s equation on a rectangle with Dirichlet boundary 
conditions, with the variables taken in the natural (typographical) ordering. 

The SOR model problem was first analyzed by Frankel in 1950 [6]. Our approach 
leads to no quantitative results that Frankel did not have, but makes it clear why the 
eigenvectors of the SOR iteration have the form they do. This analysis is restricted to 
the rectangular model problem, so it in no way supplants the much more general theory 
of matrix iterations developed by Young in the early 1950’s [10,11]. 

In 1956 Garabedian proposed a new analysis of SOR, valid in the limit of infinitesi- 
mal mesh spacing [7]. He observed that the SOR iteration is a consistent finite difference 
approximation of a time-dependent partial differential equation, so that its rate of conver- 
gence should approximate the rate of decay of solutions to that equation. To determine 
this rate, he introduced a new timelike variable s = t + x / 2 + y/2, which reduces the 
differential equation to a canonical form that can be analyzed by Fourier methods. Our 
tilting of the grid corresponds exactly to Garabedian’s introduction of the variable s. 
Thus for the SOR model problem, the consideration of a partial differential equation is 
unnecessary, and indeed the analysis in the discrete domain has the advantage that it is 
exact rather than approximate. Garabedian’s idea, however, provides additional insight 
and is applicable to more general problems. 

Approximate Fourier analysis of SOR (on the usual grid) has been discussed previ- 
ously by Kuo [8] and Chan and Elman [4] and probably others. Our tilted grid is also 
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equivalent to the “data flow times” considered by Adams and Jordan for reasons of par- 
allelizability [l]. We thank all of these authors, and also David Young, for their advice 
and encouragement. 


2. Jacobi 

Consider the discrete Poisson problem 

1 , , 

^2 ( U i-M + U i+M + u i,fc-i + u j,k + 1 “ 4u yt) = fjk, 

u jk = Fjt, 

on the square [0 ,tt] 2 , with h = ir/N. Let denote the approximation to the exact 
discrete solution u of (1) at the nth step of an iteration, with corresponding error = 
“i* ~ u ik- 

The Jacobi iteration is an example in which Fourier analysis works straightforwardly 
[ 10 , 11 ]. The errors evolve according to 

v ”* +1 = \( v j-i,k + v j+i,k + v j,k- 1 + v j,k+ 1 )» l<j,k<N-l, 

( 2 ) 

v? k =0, j = 0 ,N or k = 0, N. 

Let us consider what solutions of the form -Ik = 0(£,T7) n e t ^ x+T?v ) this iteration admits if 
we ignore the boundary conditions. We obtain immediately 

g{t,ri) = + e^ h + + e ir)h ) = §(cos £h + cos t)h). ( 3 ) 

This is the amplification factor function for the Jacobi iteration. The essential property 
is that it is an even function of £ and t)\ 

= </(-£,*?) = 9 {Z>-V) = 9 {-Z>-*l)- ( 4 ) 

If we take as initial data the linear combination 

sinfx sinryy = - c *(-«*+w) - (5) 

where £ and rj are integers in the range 1 < £ y ri < N — 1 > then the homogeneous boundary 
conditions are satisfied at n = 0 . By ( 4 ), it follows that g(£>r)) n sin£x sinrjy satisfies both 
the interior formula and the boundary conditions for all n > 0, and therefore sinfzsinrjy 
is an eigenvector of the Jacobi iteration with eigenvalue £(£,17). Since there are ( N — l) 2 


1 <j,k < N - 1 , 
j = 0, N or k = 0, N 
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of these functions, and they are linearly independent, they consitute a basis for the set 
of all grid functions {vy*}. Therefore the asymptotic convergence factor for the Jacobi 
iteration is exactly 

Jacobi _ max [I(cos£h + COSf?h)|. 

The maximum is attained with f , tj = ±1 or f , r) = ±(N — 1), 

p J ac°bi = cosh w 1 _ 1^2. (6) 


3. SOR and Gauss-Seidel 

If we attempt the same analysis for the SOR iteration, 

<#’ = (i - »K* + 

the result is 


( 7 ) 


= (!-")+ + <=—"■) + + e‘" h ), 

that is, 

l-w + + 

Now, (4) no longer holds. Therefore <jr(£,rj) does not give us the eigenvalues of the SOR 
iteration. If we find f and rj to maximize |<?(£, rj)[, there is no reason to expect the resulting 
number to describe the convergence of SOR. As it turns out, this approach produces the 
correct optimal w, to leading order in h, but a convergence rate that is too pessimistic by 
a factor of four [4]. 

Figure 1 shows how the situation can be rescued. Think of the SOR iteration as 
inhabiting a regular grid in two space and one time dimensions ( j,k,n ). Its stencil con- 
nects six points in an asymmetrical pattern, or four points in the one-dimensional case 
portrayed in the figure. Because of this asymmetry, (4) does not hold. But if we introduce 
the new “time” index 

v = 2n + j + k, (9) 

the stencil becomes symmetrical. Let us look for solutions to (7) of the form 

<fc=y(e.»?)V ( «* + ^. (10) 
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Figure 1. The SOR stencil superimposed on a space-time grid (one space di- 
mension). Introduction of the “tilted” time index v makes the stencil symmetric, 
so that Fourier analysis can be applied. The red- black labels are explained in 
Section 4. 


In the j,k y i/ variables (7) becomes 


v jk - (1 - u)v jk + -(v,._ M + v j+l k + v j k _ t + v j k+l ), 
and so a suitable value for g(£ y rj) is either root of the quadratic equation 


(ii) 


$(£> r ?) 2 = (1 - «) + ^(cos$h + cos Qh)g(S, ri). (12) 

For each £ and rj we now have a pair of amplification factors g±{£,*l), and they satisfy 
the symmetry condition (4). Therefore for any integers £, q in the range 1 < £ } r} < N — l, 
the functions g(£, T 7 )‘'sin£xsin? 7 y are eigenmodes of the SOR iteration in the v direction. 
To speak terms of eigenvectors, we note that SOR is a two-step formula with respect to 
v , but it can be recast as a one-step iteration v*' -1 ) (v*', v l/+1 ) with eigenvectors 

(sinfxsinrjy , y(£,» 7 )sin£xsin» 7 y) and eigenvalues <7(f,»?) 2 . 

It takes two steps in 1 / to advance one step in n. We conclude that the asymptotic 
convergence factor for SOR is exactly 

/>E° R = 1< ^ox max|y±(€,»?)| 2 . (13) 


In the original j,k,n coordinates, the eigenmodes become 


= y(£,r 7 ) 2n+,+k sin£xsin» 7 y, 
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and the corresponding eigenvectors are 


y (£, q) i+t sin£x sinr/y . 


(14) 


This matches the results of Frankel and others derived by different means. 

All that remains is algebra. For any £, rj, and oj > 1, the solutions to (12) are 

9±(£, *?) = ± yja 2 - (oj - 1), a = ^(cos £h + cos r]h), (15) 

and the larger of these two numbers has magnitude 


max |y± (£,»/) | 

I » 


|oc| + \Jot 2 — (a> — l) if a 2 > oj — 1, 

y/uj — 1 if O? < OJ — 1 . 


(16) 


For fixed w, this quantity can evidently be maximized with respect to £ and rj by taking 
£ = rj = 1 (among other values) and a = |cos h. 

The Gauss-Seidel iteration corresponds to w = 1. In this case (16) becomes 2|a| = 
cos h, so by (13) we have 

p GS = cos 2 h pa 1 - h, 2 . (17) 


To find the optimal overrelaxation factor for SOR, we examine the dependence of (16) 
on oj with £ = q = 1. It is readily verified that values oj > 2 lead to p^ OR > 1, so they are 
out of the running, and we can assume 1 < w < 2. In this range the second line of (16) 
obviously increases with oj , and differentiation confirms that the first line decreases with 
oj. Therefore the optimal oj is the crossover value oj — 1 = a 2 = (^cos/i) 2 , which reduces 
after a little work to 

*"-“r ^“ 2 - 2k - (18) 

By (13) and (16), the corresponding convergence rate is 


.SOR _ 1 _ 1 ~ 3m h , ol 

Popt -“opt 1 - ! + sinh « 2h - 


(19) 
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4. Relating various methods 


The change to tilted coordinates has the additional advantage of clarifying the re- 
lationships between convergence rates of different iterative methods. For example, it is 
well known that the Gauss-Seidel iteration is twice as fast as Jacobi, as is confirmed by 
comparing (6) and (17). The tilted coordinates provide a simple explanation of why this 
is so. Gauss-Seidel corresponds to the case u> = 1 of (11), and this is precisely the Jacobi 
iteration with respect to v. The factor of two comes from the fact that it takes two steps 
in v to advance one step in n . 

As another example, consider the SOR iteration with the grid points taken in the 
red-black or checkerboard ordering. This means that the Vjk with j + k even (red points) 
are processed before the t;y* with j + k odd (black points), and the iteration takes the 
form 


t # 1 = (i-"K* + 

v y * +1 = (1 ~ + 


1 Wj-l,k + v ?+l,k + v ?,k-l + v ?,k+l) 


for j + k even, 
for j + k odd. 


For each a ;, this method has the same convergence rate as the iteration (7) in the natural 
ordering. Young proved this algebraically by determining the eigenvalues and eigenvectors 
of the associated iteration matrices [11]. Again, the change to tilted coordinates gives a 
more intuitive explanation, as illustrated in Figure 1 for the case of one space dimension. 
At step v we are computing only the red points, say, and at step v+1 only the black points. 
Recasting SOR as a one-step iteration (t;*^ 2 , v 1 ' -1 ) ^ (v”, v 1 '"*’ 1 ), as in the last section, 
we obtain simply the red-black ordering. Thus Figure 1 can be viewed as depicting an 
SOR iteration either in j 9 n coordinates with the natural ordering, or in j y v coordinates 
with the red-black ordering. Hence these two orderings must have the same asymptotic 
convergence rate. 


The conclusions above depend on the fact that the convergence rate is independent 
of the particular initial data used, depending only on the eigenvectors. Note that we can 
switch back and forth between arbitrary data at fixed n or at fixed v by taking partial 
iterations over a triangular portion of the grid. In fact, writing out these partial iterations 
algebraically gives a similarity transformation relating the iteration matrices. 

In this paper we have considered just the five-point Poisson model problem, and 
presented an easy way to obtain classical results with, we hope, additional insight. The 
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tilted grid may also prove useful in obtaining new results. It has already been applied to 
settle a conjecture of Adams and Jordan regarding the equivalence of certain orderings 
for the nine-point stencil [1]. These results will be reported in [2]. 


References 

[1] L. M. Adams and H. F. Jordan, Is SOR color-blind?, SIAM J. Sci. Stat. Comp. 7 
(1986), 490-506. 

[2] L. M. Adams, R. J. LeVeque, and D. M. Young, Analysis of the SOR iteration for the 
9-point Laplacian, to appear as an ICASE report. 

[3] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comp. 
31 (1977), 333-390. 

[4] T. F. Chan and H. C. Elman, Fourier analysis of preconditioned iterative methods, 
Yale Univ. Dept, of Comp. Sci. Rep. Yaleu/DCS/RR-410, to appear. 

[5] G. E. Forsythe and W. R. Wasow, Finite- Difference Methods for Partial Differential 
Equations, Wiley, 1960. 

[6] S. Frankel, Convergence rates of iterative treatments of partial differential equations, 
Math. Comp. 4 (1950), 65-75. 

[7] P. Garabedian, Estimation of the relaxation factor for small mesh size, Math. Comp. 

10 (1956), 183-185. 

[8] C.-C. Kuo, B. C. Levy, and B. R. Musicus, A local relaxation method for solving 
elliptic PDEs on mesh-connected processor arrays, SIAM J. Sci. Stat. Comp., to 
appear. 

[9] R. D. Richtmyer and K. W. Morton, Difference Methods for Initial- Value Problems, 
2nd ed., Wiley-Interscience, 1967. 

[10] R. S. Varga, Matrix Iterative Analysis, Prentice-Hall, 1962. 

[11] D. Young, Iterative Solution of Large Linear Systems, Academic Press, 1971. 


7 



Standard Bibliographic Page 


1. Report No. NASA CR— 178191 2. Government Accession No. 

ICASE Report No. 86-63 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

FOURIER ANALYSIS OF THE SOR ITERATION 

5. Report Date 

September I 986 

6. Performing Organization Code 

7. Author(s) 

Randall J. LeVeque, Lloyd N. Trefethen 

8. Performing Organization Report No. 

86-63 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665-5225 

10. Work Unit No. 

11. Contract or Grant No. 

NAS1-18107 

13. Type of Report and Period Covered 

C r\r\ frnpfnr Ponnrf- 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

cnr r\r\ m ai 

15. Supplementary Notes J ° Ui 


Langley Technical Monitor: SIAM Num. Analysis 

J. C. South 


Fi nal Report 

16. Abstract 


The SOR iteration for solving linear systems of equations depends upon an 
overrelaxation factor 1 . We show that for the standard model problem of 
Poisson's equation on a rectangle, the optimal 1 and corresponding convergence 
rate can be rigorously obtained by Fourier analysis. The trick is to tilt the 
space-time grid so that the SOR stencil becomes symmetrical. The tilted grid 
also gives insight into the relation between convergence rates of several 
variants. 


17. Key Words (Suggested by Authors(s)) 


18. Distribution Statement 


successive over relaxation, SOR, 


64 - Numerical Analysis 


iterative methods, Fourier analysis 

67 - Theoretical Mathematics 



Unclassif iec 

1 - unlimited 

19. Security Classif.(of this report) 

20. Security Classif.(of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 

9 

A02 


For sale by the National Technical Information Service, Springfield, Virginia 22161 
NASA Langley Form 63 (June 1985) 


NASA-Langley, 1986 




